#! /usr/bin/env python
"""
median editor

Copyright (C) 2012-2013 Christian T. Steigies <steigies@physik.uni-kiel.de>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""

__author__ = "Christian T. Steigies <steigies@physik.uni-kiel.de>"
__license__ = "GPL License"

## --------------------------------------------------------------------------
def median(k, n0, e0, nt):
    """
    k      - number of channels
    n0[k]  - initial count rates (baseline)
    e0[k]  - initial efficiencies
    nt[k]  - current countrates 
    ne[k]  - corrected countrates
    et[k]  - calculated efficiency
    """
    
    i = j = ik = jk = ii = 0
    c1 = c2 = s0 = sn = tet = 0.0
    # define all lists seperately, otherwise the same list will have different names!
    ne = [0.]*k
    et = [0.]*k
    r = [0.]*k
    
    for i in range(0, k):
        if nt[i] > 0.0:
            ik += 1 # number of counters with countrate >0
        s0 += n0[i]/e0[i] # sum, corrected by efficiency
        r[i] = nt[i] / e0[i] / n0[i] # ratio countrate versus reference value, corr by eff
        et[i] = r[i] # save efficiency for each counter

    if ik < 3: # if less than three counters, return 0, EXIT
        median_sum = 0.0
        for i in range(0, k):
            et[i] = 0.0

    else: # three or more counters
        jk = ik/2 + 1 # Note: Ik DIV 2 + 1
        ## find median of ratio
        for j in range(0, jk):
            c1 = 0.0
            for i in range(0, k): # find max ratio and index
                if r[i] >= c1:
                    c1 = r[i] 
                    ii = i
            if (j+1) < jk: # for last value (jk) set ratio of max to 0.0
                c2 = c1
                r[ii] = 0.0 

        tet = c1 # ratio of max
        jk = (ik / 2) *2    
        if jk == ik: # even number of counters with countrate > 0.0
            tet = (c1+c2)/2.0 # tet = average of ratio of max and second last counter
    
        median_sum = 0.0
        for i in range(0, k):
            et[i] /= tet
            if et[i] != 0:
                ne[i] = nt[i] / et[i]
            ne[i] = int(ne[i] + 0.5) # return counts as integers
            median_sum += ne[i] # or medianSum = s0*tet

    return(int(median_sum+0.5), ne, et) # sum, channels, efficiency
